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ABSTRACT 

A simple rapid method is described for determining the spectrum of a surface field 
(in spherical harmonics) from harmonic analysis of direct (in situ) measurements along great 
circle arcs. The method is shown to give excellent overall trenu r (smoothed spectra) to very 
high degree from even a few short arcs of satellite data. Three examples are taken with 
perfect measurements of satellite tracking over a planet made up of hundreds of point-masses 
using (1) altimetric heights from a low orbiting spacecraft. (2) velocity (range rate) residuals 
between a low and a high satellite in circular orbits, and (3) range-rate data between a station 
at infinity and a satellite in a highly eccentric orbit. 

In particular, the smoothed spectrum of the Earth’s gravitational field is determined to 
about degree 400(50 km half wavelength) from i° X 1° gravimetry and the equivalent of 
1 1 revolutions of Geos 3 and Skylab altimetry. This measurement shows there is about 
46 cm of geoid height (rms world-wide) remaining in the field beyond degree 1 80. 


•On leave of absence from the Universit) of New South Wales, Sydney, Australia. 
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GRAVITATIONAL SPECTRA FROM DIRECT MEASUREMENTS 


INTRODUCTION 

The spectrum of a two-dimensional field in spherical harmonics is conveniently ex- 
pressed by the so-called degree variances of its fluctuations. We demonstrate here that the 
trend of the spectrum is far easier to find than the actual field. For example the anomalous 
gravitational potential external to a planet can be written as: 


¥■ 2 m-0 ' ' 


p lim (sin 


Cg cos mX + S^ m sin mX 


(I) 


where r is the distance to the center of mass. R is the radius of the smallest sphere enclosing 
all the mass. 0 and X are planet centered latitude and longitude, v is the planet’s gaussian 
mass, P^ m are fully normalized associated legendre polynomials (lleiskanen and Moritz. I%7. 
p. 31) and C^ m . S» are their associated coefficients. The potential degree variances are: 


m - 0 ' ' 


O 


We note that this power spectrum (and others simply related to it) is the same for any 
choice of reference pole for the field. Potential fluctuations on any surface great circle 
always has maximum wave number C with respect to the ( 2d + I ) terms of that ‘line’ of the 
spectrum. 

I he internal density source lot both the external field <>i its Minimal y s|m i tiiiin is »d 
course not unique. But the mass integral representation of the field coefficients together 
with other physical evidence puts important constraints on the source spectrum. From 
llotinc, 1 9b9, p. I 5*E 
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d\, (3) 


whore d A is an a ca element, p is the density and p is the Earth’s average density. For 
example, the severe radius factor in (3) means that the deep density anomalies can have little 
affect on the high degree field if they are of reasonable sire (e.g., Jeffreys, 1962) or are not 
to give rise to unreasonably large shearing stresses (e.g., Higbie and Stacey, 1971a, b). 

Kaula ( 1966, p. 9«) first gave the simplest expression of the Earth’s gravitational spec- 
trum ( from autocovariance analysis of gravimetry): 


o“ = (2C+ l)(10 - 5 /«V. 

a result which we will later show is a fair overall approximation, though too high up to 
about degree 60 

in recent years many attempts have been made to interpret the source of this spectrum 
using mineralogic. scismologic, thermal as well as stress evidence (e.g., Kaula, 1963; Higbie 
and Stacey. 1971b; Bott. 1971; Allan, 1972; La ni beck. 1976; Dzicwonski, 1977). Kaula, 
1977 concludes that the smoothness of the spectrum is probably as significant as its general 
form (This smoothness may be due to ‘dynamic’ compensation at depth as a result of 
convective processes.) 

Our goal here is more modest; to show that simple harmonic analysis of ‘central’ arcs 
of ‘field’ information will rapidly elucidate the spectrum. Our method complements Kaula’s 
( 1959) classical technique of autocovariant analysis (see also Hciskanen and Moritz, 1967, 
pp. 252-259; r scheming and Rapp, 1 974). But the new method does not require global 
information and is more rapid (with access to the fast Fourier transform). 
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METHOD 


We observe from (1) that harmonic analysis of potential data on the equator or any 
latitude circle of a spherical surface will yield Fourier coefficients [for wave numbers n * m] 



C £n* S £n 


(4) 


The expected value of the one dimensional power spectrum for this circle data is 


E(j2 +T 2 ) * (- 

v n,co* n, sin yr 



R \ 2 ' S 2 


P en (sin*) 


t»Jl 


EC jL + EE en 


(5) 


assuming that the field coefficients are chosen from an uncorrelated population. But this is 
a weak assumption and in any case we really seek the average power spectrum for a given 
field over all possible great circles (or small circles of the same latitude). Consider the case 
for great circle data. We assert (without proof here) that: 
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v n.cos n.sin' 


( 6 ) 


where 


2 x 
* ' (28 + 1) 


the mean square coefficient by degree. 

Note that for coefficients chosen from an uncorrelated population with 
ECg m = , (6) follows from (5) because the degree variance of an harmonic field is 

invariant under rotation. We choose the equator to evaluate the associated legendre 
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polynomials because of the simplicity and stability of the formulas for them there. For 
example, from Abramowitz and Stegun, 1964, we find: 


S + 2,n 


(o) - 


(8~n + l)(fi + n+ l)(2fi + 5) 


(fi - n + 2) (8 + n + 2) 


? 8, n (0) 


( 7 ) 


Using the spherical approximation of Bruns* formula [Heiskanen and Mortiz, 1967, 
pp. 84-85] the average great circle geoid height (N) spectrum is: 


< ( N n,c M + N2 „,«n»= 2R 


—(fi-n even) 
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( 8 ) 


8 - n 


Since the gravity anomaly field at the surface of the Earth is given by [e.g., Heiskanen and 
Moritz, 1967, p. 89] : 
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(8- l) P gm (sin £) 


$-2 ni = 0 


C" cos mX + S„ sin mX ; 
8rn 8m 


the average of the one dimensional power spectra for all great circle arcs of gravity anomaly 
data is. 


• (6 -neven) 
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n.cos n.sin 
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(9) 


where 



Similarly it is easily shown {from the radial derivative of (1) and comparison with (6)] 
that the average spectrum for anomalous fadial acceleration external to the sphere is: 
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( 10 ) 
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Of more general interest (than the radial) is (he anomalous line-of-sight acceleration 
with respect to a fixed external direction. For example the (t?$i Jual) tracking data returned 
from space probes or through high-low satellite-satellite exchange contains this kind of 
information. The line-of-sight of the signal actually changes slowly during the tracking of 
the probe or low satellite due to the planet’s rotation or the fact that the high satellite is not 
at infinity. But these defects actually enhance tne averaging of aspects necessary to achieve 
good global results. 

The gravitational acceleration in the direction of the planet's spin axis (Z) provides a 
convenient starting point for the “fixed” direction spectrum. From Hotine (1969. p. 180, 
21.140/141) we find these accelerations are: 
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where 
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( 12 ) 


Again, the average of all great-circle spectra of this (single direction) field can be expressed 
simply in terms of the mean square coefficient by degree for (a, b) as in (6). (9) and (10). 
Furthermore, we can show from ( 1 2) that this variance is: 


6J(a.b) 



(13) 



Situs Sq _ j is invariant for any position of the z axis, the average of every fixed direc- 
tion <D) spectrum is the same. We find: 

••(f -n even) 

<(D 2 + D 2 . ) > = 2y 2 V* 
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(14) 


Comparing (14) to (10) the average fixed direction spectrum of accelerations is no less than 
2/3 of the radial for realistic smooth gravitational spectra such as Kaula's rule. However, if 
the probe is observed entirely cross track (normal to its orbit plane) the acceleration spec- 
trum is weaker still since the leading tesseral harmonic is absent. In this case only, ( 12) 
gives the average cross track (CT) spectrum directly (for all “equatorial" great circle tracks 
observe ! from a pole; 0 = 0): 
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Finally, the data type most frequently encountered in spaa* tracking is range rate or 
relative velocity between the observer and the probe. Since acceleration (A) is the derivative 
of velocity (V) we have used tne harmonics of tracks of V data to estimate the corresponding 


wave numbers of the A spectrum; 


| A - + A‘ • 

1 n.cos n.sm 



V- 

n.cos 




V : j 

n.sm 1 


( 10 ) 


where P is the period of the probe's orbit. 
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Inversion 


Imagine (hat the average circular “track" power spectra (Pj^ ) are available for a certain 
data type. All of these are related to an infinite sum of every other degree variance of the 
spherical harmonic field starting at the degree equal to the wave number of the track spectra 
(or one less for “fixed" direction spectra). 

•*(e-n even) 

P„ 2 = ^ S 2 (8, n) |Sj or 8j,,l (17) 

e - n 

where S" are sensitivity factors for the various data types calculable from (8), (9), (10), (14) 
or (16). 

In practice this infinite system is solved iteratively for the spherical harmonic variances 
5g from high to low degree. It starts from a sufficiently high degree beyond the significant 
power spectrum so that the truncation error (from assuming 6^ ^ = 0) is negligible. 

In the “real world” probes do not follow great circles over planets nor do observers 
view them from fixed directions over complete (global arcs). But global power from .ion- 
global arcs can be estimated from well known data “windows” Furthermore, trajectories 
of close satellites of even rapidly rotating planets are often close enough to great circles and 
observing directions constant enough for this simple methodology to be useful. 
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SIMULATIONS 


We have constructed a model (spherical) Earth, rotating with 266 randomly chosen 
point masses added in three layers (at 20.1 50 and 1000 km depth) so as to give fairly realistic 
degree variances (figure 1 ). In the first simulation, ten (1/6 revolution) arcs of perfect gcoid 
height altimetric data (every 10 seconds) from a 200 km circular orbiting satellite (115° 
inclination) has been taken (figure 2). The average track-power spectrum (rins) for these 
arcs are given in figure 3 together with the values expected from the degree variances for 
this model [calculated from (3)1 using (o). The fidelity of these averages is remarkable 
considering the field’s high local roughness and its limited sampling here. The global 
power estimates were made from residuals to simple end-matching trend lines, using a transfer 
function P~(global) = P“( measured )/6, for the 1/6 revolution data | Wagner. 1977, Appendix 
C|. A smooth power law of the form An® * ( * lo>: n * has been fitted to these averages and 
smooth mean square coefficients estimated from back-solving equation (8). The results 
( figure I ) confirm the overall accuracy of the method even to very high degree. 

Spectral analysis has been made of four ( 1/2 revolution) arcs of perfect range rate data 
(taken every 10 seconds) between a st. tionary satellite at 6.61 Earth radii and the 200 km 
satellite (figure 2). Average velocity power at altitude falls off steeply (figure 4). Starting 
at about 60 cycles, the leakage or edge effect of the “rectangular window” begins to domi- 
nate the spectrum. Using residuals to quadratic trend curves which match the first derivatives 
as well as the data at the ends, the average power spectrum is dramatically improved to about 
170 cycles. Some aliasing and leakage is evident at still higher frequencies. 

Again a smooth power function has been lit to the average spectral data (using the 
linear trend results up to 60 cycles and the quadratic to 166 cycles). Inversion of this func- 
tion | from (14) and (16)1 is reasonably accurate up to about 100 cycles (figure I ). Above 100 


cycles the smooth variances art systematically high, « result due to the lack of coverage of 
the back side of the planet. 

It seems likely that unless accuracies within 10% are required, it will not be necessary 
to correct for Earth rotation or the 1 T variation of direction to the low satellite from 
synchronous orbit. 

Finally, we have a promising result using a single pass of simulated range rate data from 
a Venus-orbiter type mission. Here we tracked a highly eccentric 24-hour orbiter from a 
truly fixed direction normal to the orbit plane. We used the same model planet as before 
but without rotation (figure 2). The significant data is highly compressed (in time) around 
periapsis (200 km altitude) and our analysis is of a stretched out version of the data at equal 
mtervak of true anomaly (central angle & ). Figure 5 shows the line-of-sight acceleration of 
this probe between ± 30° of periapsis (data every .36® ). The power spectrum (using the 
“rectangular” window) shows high frequency leakage beginning at about 1 50 cycles, but is 
otherwise quite regular ( figure 6). We have estimated the spectrum for a circular orbit at 
various radii along this trajectory [using (14) and (1 5)1 but no single spectrum or average of 
spectra “fits” the trend of this data satisfactorily over all frequencies. However, the expected 
circular orbit cross track spectrum is an improvement over the corresponding general fixed 
direction spectrum at high frequencies ( n > 50). Examination of trajectories at different 
longitudes shows that while the mid-frequency spectrum in figure 6 is indeed higher than 
average, no simple circular orbit spectrum can accurately represent all frequencies. Neverthe- 
less, the circular orbit spectra are fairly good first approximations. 
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SMOOTH SPECTRUM OF THE GEOPOTENTIAL 

We have determined the smooth spectrum of the geopotential to high degree from 
average track spectra using: (1) satellite altimetry over the oceans and (2) 1° X 1° global 
gravimetry. There is remarkable agreement between these two independent methods. 

The altimetry consists of 8 1 arcs of sea surface heights (R. Hanley, private communi- 
cations, 1077; McGoogan, Leitao and Wells, 1975) over all the world’s oceans estimated by 
the radar altimeter aboard the Skylab IV and Geos 3 spacecrafts (figure 7). The individual 
heights represent averages over 1 5 to 20 km but a few arcs have teen sampled at up to 70 km 
spacing (Wagner, 1977; Wagner, 1978). The original data is noisy (±15 to 50 cm) and is 
contamina with low frequency orbit error (1 to 10 m; n “ 0, l, 2 cycles/rev.). The orbit 
error and most of the geoid signal to 20 cycles has been removed by fitting the heights to the 
geoid from Goddard Earth Model 7 (Wagner et al,. 1977) and empirical terms with these low 
frequencies. Group averages of the (Gem 7) residual spectra for these arcs are shown in 
figure 8, The “noise" tail starting at about 750 cycles corresponds to an average “white 
noise” level of 2 1 cm. 

Wc have also evaluated the spectrum of sea surface topography (departures from the 
level geoid). Topography which changes with time has been evaluated from ( 1): eleven 
Geos 3 altimeter arc overlaps (within 10 km. having I *4 and 214 month time separation); 

( 2 ) a number ol profiles over sea level air pressure maps (c.g., Straliler, l%7) interpreted as 
an inverted barometer on the sea surface; (3) a 20-mmute Geos 3 profile in the Pacific of 
luni-solur tide heights ( Estes, 1977); and (4) similar profiles over bi-weekly averages of 
topography for the East Australian current (Boland and Hamon, 1 970). Actually (2) and (4) 
include “permanent" effects ir. the lower frequencies. 

The "permanent" topography has been estimated from (1), four long profiles of Defant’s 
( 1 9(> 1 ) map of the North Atlantic (in; Apel, 1976} and (2). five similar profiles ever annual 
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averages of surface topography for the Pacific (Wyrtki, 1975). Croup averages for these 
spectra are displayed on figure 8. Note that at frequencies higher than 50 cycles the dynamic 
(time varying) topography dominates the departure spectrum. The altimeter-measured 
departure spectrum is fairly flat from 50 to 500 cycles. But above 500 cycles we have not 
evaluated enough data to properly distinguish it from the (noise-reduced) surface height 
spectrum. Using a composite average of these departure spectra and subtracting from the 
noise-reduced sea height spectrum, an estimated smooth geoid-track spectrum has been 
derived. 

The inversion of the smooth geoid-track spectrum according to (8) is shown in figure 9 
along with the same result from Wagner. 1978 using only about 1/2 of this data. The 
stability of the solution appears to be excellent. It is also confirmed by separate solutions 
from independent data-arc sets. 

An entirely different data set. 58,000 ! n X 1° mean gravity anomalies (R. II. Rapp, 
private communication. 1977) has also been analyzed by this method. The remaining 26.800 
global 1° X 1° anomalies have been estimated by averaging over the 8 adjacent 1° blocks 
centered on the missing data. All missing blocks were first filled in with 5° X 5" mean 
values (Rapp. 1977). 

Harmonic analysis of the 180 global meridional arcs of this data set was performed with 
the fast I'ourier transform (Oppenheim and Schafer. 1975). The average power spectrum 
from this analysis is shown in figure 10. 1 he overall comparison is excellent with the spec- 
trum expected from the previous analysis of marine altimetry. In making this comparison it 
is essential to use ( in (9)| geopotential variances reduced for the effects of 1° X 1° smoothing 
(Pcllinen. 1966; Rapp. 1 977. pp. 18-20). The anomaly power remains high for high degree, 
and the track spectra at n are sensitive to all geopotential variances above n of the same 
parity. Thus the smoothing effect on 1° anomalies is seen to extend back as far as 10 cycles. 


NSe interpret the excess of track power above 100 cycles to be due principally to aliasing 
because no “noise tail” is obvious. (It should be present at an average white noise level above 
12 mgals.) 

NVe estimate tlu aliasing in a declining power spectrum, at and below the Nyquist limit 
(here. N ISO cycles) to be given by: 


(P*).. , 

n aliased 


( I* ) 

ii expected 


I + 


(~r*(v-f 


(IS) 


where, 


( P- ) , . 

u expected 


(An“ B ) : , B >1/2 . 


(l l >) 


I'he lower limit on the slope of the “true” spectrum (B) ensures that the total signal is 
finite. 

I rom figure 10, we estimate that B i is the correct slope of the spectrum near 
the Nyquist limit, somewhat flatter above and steeper below. ( 1 he final result is fairly 
insensitive to this estimate. N)e also note that for this slope the effect of the second fold 
of frequencies is only I 1 of the first at n N.) Using (IS) we have corrected the 
average track spectrum for aliasing. I hen. using the odd and even frequencies separately 
(in conjunction with a smooth extrapolation of the power beyond n ISO to account for 
the truncation) we have found both odd and even degree variances from inversion of ( l >). 

I he result lor individual degree variances however is too widely vary ing to be realistic 
Recall that the track spectrum should be smoother than the field spectrum since it is formal 
Iront in infinite sum ot the degree variances. Evidently the meridional arcs do not yield a 
tan sampling ol the field at all frequencies, only showing north-south variations anil over- 
emphasizing the polar regions. Nevertheless, averaging the variances in It) degree groups 
( t ■> 20) gives a smoot lie i result ( figure ‘>1 agreeing remarkably well with the marine geoid 
spectrum ( V ^ 25) and ( V! < 20) the recent tioddard Earth Model 1 0B il erch et al., U)78). 
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Also shown on figure 9 is an early mollification of Kaula’s rule from covariance analysis 
of gravimetry (Kaula. I%9). The significant change from 10 ~ 5 /8‘ occurs for 8 >90 showing 
greater power than the “rule,” which we confirm from 8 > 60 to at least 8 = 400. 

In summary, the gravitational spectrum of the Harth | from extensive harmonic analysis 
of track spectra of marine altimetry (If > 20). surface gravimetry (8 < ISO) as well as from 
least squares solutions for satellite tracking, anil broader altimetry and gravimetry data 
(8 < 3b)| is given closely by: 


8 , = 


where: 


10 6 a= 12.6 . 

b* 2.198 ; 

8 < 10 

10 6 a«50.1 . 

b = 2.797 ; 

10 < 6 < 20 

I0 6 a = 0.430, 

b = 1.209 ; 

20 < 8 < 50 

10^1 = 2.74 . 

b = 1.681 ; 

50 < 8 <. 1 50 

I0 h a = 18.5 . 

b = 2.063 ; 

1 50 < 8 < 400 

O' 

u 

II 

t J 

b = 2.1 ; 

8 > 400 


This spectrum has been used to calculate the total geoid power above each complete 
degree (ruts, world-wide): 


IV I B 1 

(. max 


The same power from Kaula’s rule is closely approximated by (64^ ) m 

for 8 mtx > 10 (e g.. Chovit/.. ! l >73). Figure I I shows that even if the geopotential is known 
perfectly to (ISO. ISO) about do cm ot power remains in the geoid undulations. It should be 
emphasized that the piece w is; spectra of (21) are smooth and that variations for individual 
degrees can be as great as 2S‘< (see figure 9). 
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SUMMARY AND CONCLUSIONS 

An efficient method has been found to determine spherical harmonic spectra directly 
from harmonic analysis of track data. The method is a ‘‘short cut” complement to the com* 
plete covariance analysis of data on a sphere. It does not require global data and can make 
use of the fast Fourier transform. The method has been developed for a wide variety of 
satellite tracking (as well as surface) data types for application to the determination of 
gravitational spectra. Included is a simple algorithm to allow for aliasing. In particular, the 
smooth spectrum of the geopotential has been derived from global surface gravimetry and 
the equivalent of 1 1 revolutions of marine satellite altimery. Important points to note are: 

( 1 ) The present application is only strictly for probe data on great circle tracks; 

(2) if the tracks are global, individual degree variances can be found but their accuracy 
depends on sufficient number and distribution of tracks; 

(3) if the tracks are not global, smooth (degree-averaged) spectra are found; 

(4) results from non-global tracks within about 10% of great circles show smooth 
spectra with no significant distortions; 

(5) results from highly elliptical tracks (rectified) do show significant distortions but 
may still be useful. 
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Figure 1. RMS Gravitational Coefficient by Degree, for Model Earth 
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Figure 3. Average Power Spectrum for 10 Arcs of Altimetry over Model Earth 
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Figure 4. Average Power (rms) in Range Rate between High and Low Orbite 

over Model Earth 
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Figure 5. Acceleration of the Line of Sight between an Eccentric Orbitcr and 

a ‘Point at Infinity’ 
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Figure 7. Distribution of Altimetry (HI Arcs and 1 1 Overlap Pairs) in Evaluation of Geoid Spectrum 
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Figure X, Sea Surface Power Spectrum (nils) from Satellite Altimetry and 

Oceanographic Data* 
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Figure 9. Geopotential Spectrum* 
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